#include <stdio.h>
#include <math.h>

typedef	float VECTOR[3];
VECTOR	vel[600][248][248];
float	curl_mag[600][248][248];

main(unsigned argc, const char *argv[])
{
	if (argc != 2) {
		fprintf(stderr, "Usage: %s INPUT_FILE\n", argv[0]);
		fprintf(stderr, "      INPUT_FILE: For example, velocity.0030.txt\n");
		return -1;
	}

	// Figure out what input file to open and open it.
	const char *input_file_name = argv[1];
	FILE *infile;
	if ( (infile = fopen(input_file_name, "r")) == NULL) {
		fprintf(stderr,"Could not open %s for reading\n", input_file_name);
		return -2;
	}

	// Figure out what VTK output file to open and open it.
	char	v_output_file_name[4096];
	sprintf(v_output_file_name, "%s.curlmag.vtk", input_file_name);
	FILE *v_outfile;
	if ( (v_outfile = fopen(v_output_file_name, "w")) == NULL) {
		fprintf(stderr,"Could not open %s for writing\n", v_output_file_name);
		return -3;
	}


	// Figure out what binary output file to open and open it.
	char	b_output_file_name[4096];
	sprintf(b_output_file_name, "%s.curlmag.bin", input_file_name);
	FILE *b_outfile;
	if ( (b_outfile = fopen(b_output_file_name, "wb")) == NULL) {
		fprintf(stderr,"Could not open %s for writing\n", b_output_file_name);
		return -3;
	}

	// Scan the input file into the velocity field.
	printf("Reading velocity file\n");
	unsigned x, y, z;
	for (z = 0; z < 248; z++) {
	  for (y = 0; y < 248; y++) {
	    for (x = 0; x < 600; x++) {
		if (fscanf(infile, "%f %f %f", &vel[x][y][z][0], &vel[x][y][z][1], &vel[x][y][z][2]) != 3) {
			fprintf(stderr,"Could not read cell (%d,%d,%d)\n", x,y,z);
			return -4;
		}
	    }
	  }
	}

	// Compute the curl field.
	printf("Computing Curl\n");
	for (z = 0; z < 248; z++) {
	  for (y = 0; y < 248; y++) {
	    for (x = 0; x < 600; x++) {
		// Zero the entries at the upper end of
		// the X, Y and Z dimensions; there is no way to compute them (no upper
		// cell value to use).
		if ( (x == 599) || (y == 247) || (z == 247) ) {
			curl_mag[x][y][z] = 0.0;
		} else {
			VECTOR	curl;

			curl[0] = ( vel[x][y+1][z][2] - vel[x][y][z][2] - vel[x][y][z+1][1] + vel[x][y][z][1] ) / 0.001;
			curl[1] = ( vel[x][y][z+1][0] - vel[x][y][z][0] - vel[x+1][y][z][2] + vel[x][y][z][2] ) / 0.001;
			curl[2] = ( vel[x+1][y][z][1] - vel[x][y][z][1] - vel[x][y+1][z][0] + vel[x][y][z][0] ) / 0.001;

			curl_mag[x][y][z] = sqrt(curl[0]*curl[0] + curl[1]*curl[1] + curl[2]*curl[2]);

		}

	    }
	  }
	}

	// Write the header to the VTK output file telling what kind of VTK file
	// this is and what its data set is
	fprintf(v_outfile, "# vtk DataFile Version 2.0\n");
	fprintf(v_outfile, "Calculated Curl Magnitude\n");
	fprintf(v_outfile, "ASCII\n");
	fprintf(v_outfile, "DATASET STRUCTURED_POINTS\n");
	fprintf(v_outfile, "DIMENSIONS 600 248 248\n");
	fprintf(v_outfile, "SPACING 0.001 0.001 0.001\n");
	fprintf(v_outfile, "ORIGIN 0 0 0\n");
	fprintf(v_outfile, "POINT_DATA 36902400\n");
	fprintf(v_outfile, "SCALARS density float 1\n");
	fprintf(v_outfile, "LOOKUP_TABLE default\n");

	// Write the curl field into the VTK and binary files.
	printf("Writing VTK and binary files\n");
	for (z = 0; z < 248; z++) {
	  for (y = 0; y < 248; y++) {
	    for (x = 0; x < 600; x++) {
		if (fprintf(v_outfile, "%0.3E\n", curl_mag[x][y][z]) <= 0) {
			fprintf(stderr,"Could not write VTK cell (%d,%d,%d)\n", x,y,z);
			return -5;
		}
		if (fwrite(&curl_mag[x][y][z], sizeof(float), 1, b_outfile) != 1) {
			fprintf(stderr,"Could not write binary cell (%d,%d,%d)\n", x,y,z);
			return -6;
		}
	    }
	  }
	}

	// Done!
	fclose(infile);
	fclose(v_outfile);
	fclose(b_outfile);

	return 0;
}